Photonic analogue of Zitterbewegung in binary waveguide arrays 
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An optical analogue of Zitterbewegung (ZB), i.e. of the trembling motion of Dirac electrons caused by the 
interference between positive and negative energy states, is proposed for spatial beam propagation in binary 
waveguide arrays. In this optical system ZB is simply observable as a quiver spatial oscillatory motion of the 
beam center of mass around its mean trajectory. © 2009 Optical Society of America 
OCIS codes: 230.7370, 000.2658, 050.5298 



Zitterbewegung (ZB) refers to the rapid trembling mo- 
tion of a free Dirac electron caused by the interference 
between positive and negative energy state components 
[1,2]. However, the amplitude of ZB oscillations turns 
out to be extremely small, of the order of the Comp- 
ton wavelength, which defines the limit of electron lo- 
calization. Therefore, its direct observation is unlikely, 
though observable consequences of ZB have been found 
in the response of electrons to external fields [1]. Sim- 
ilarly to other effects, like the Klein paradox [2], ZB 
has been for long time regarded as a relativistic effect 
rooted in the Dirac equation. However, several authors 
have recently shown that ZB is not unique to Dirac elec- 
trons, rather it is a generic feature of wave packet dy- 
namics in spinor systems with certain linear dispersion 
relations, such as those exhibiting the so-called Dirac 
points (DP) that describe massless fermions [3-13]. In 
condensed-matter and matter- wave physics, trapped ions 
[3], graphene [4,5], and ultracold neutral atoms [6] have 
been proposed as candidate systems for a direct observa- 
tion of ZB. Similarly, in the optical context it was shown 
that DP can occur in two-dimensional photonic crys- 
tals [7-11] or in negative-zero-positive index metamate- 
rials [12]. Photonic analogues of ZB have been recently 
proposed in Refs. [11,13], in which the observation of ZB 
requires time-resolved measurements of pulse transmis- 
sion through photonic crystal or metamaterial slabs of 
different thickness. In this Letter we propose a simpler 
optical system, consisting of a binary waveguide array, 
which enables an easy visualization in space of photonic 
ZB. The ability of mapping typical ultrafast phenom- 
ena occurring in the matter as spatial propagation of 
light waves in coupled waveguide structures has been 
successfully demonstrated in several experiments (see, 
e.g., [14-17]), and our proposal should thus greatly fa- 
cilitate the way toward a first observation of ZB. 

Let us consider propagation of monochromatic light 
waves at wavelength A in a one-dimensional optical lat- 
tice, which is described by the following scalar equation 
for the electric field envelope E(x,z) 
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where n s is the substrate refractive index and n(x) is the 



Fig. 1. (Color online) (a) Schematic of a one-dimensional 
binary array (left) and refractive index profile (right), (b) 
Dispersion curves of the first two minibands of a tight- 
binding binary array (solid curves), and correspond- 
ing dispersion curves of the Dirac equation (4) (dotted 
curves), (c) Broad-beam excitation geometry of the ar- 
ray at a tilting angle 0. 



refractive index profile of the lattice. To realize a pho- 
tonic analogue of ZB, a two-band model with a dynamics 
described by a two-component spinor wave function is 
needed. Such a model can be realized by either a singly- 
periodic lattice with a shallow sinusoidal refractive index 
profile, where the two components of the spinor wave 
functions are related by simple linear transformation to 
the amplitudes of counterpropagating waves in the lat- 
tice (see, for instance, [18]), or to a tight-binding bi- 
nary superlattice composed by two interleaved sublat- 
tices A and B [see Fig. 1(a)], where the spinor wave func- 
tions correspond to the occupation amplitudes in the 
two sublattices. Here we study ZB in the latter struc- 
ture and assume a symmetric intersite coupling, which 
corresponds to the experimental conditions of Ref. [17]; 
the analysis could be extended, if needed, for the more 
general case of non-symmetric intersite couplings [20] . In 
the tight-binding approximation, light transport in the 
binary lattice is described by coupled-mode equations 
for the modal field amplitudes c n in the various waveg- 
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uides [17,20,21] 

i(dc n /dz) = -a(c n+ i + Cn_i) + (-l) n (5c n , (2) 

where 26 and cr are the propagation constant mismatch 
and the coupling rate between two adjacent waveguides 
of the array, respectively. The tight-binding model (2) 
supports two minibands, whose dispersion curves are 
readily obtained by the plane- wave Ansatz c n (q) ~ 
ex.p(iqna — iuoz) and read [20] 

u(q) = ±y / 5 2 + 4a 2 cos 2 (qa) (3) 

[see Fig. 1(b)]. Let us assume that the array is excited by 
a broad beam (e.g. Gaussian shaped) incident onto the 
array at an angle close to the Bragg angle = A/ (4n s a), 
i.e. E(x,0) = G(x)exp(27ri0n s x/\), where a is the spac- 
ing between adjacent waveguides and G{x) varies slowly 
on the spatial scale ~ a [see Fig. 1(c)]. At such an inci- 
dent angle, the modes in adjacent waveguides are excited 
with a nearly equal amplitude but with a phase differ- 
ence of 7r/2. After setting C2 n (z) — (— l) n, 0i(ri, z) an d 
C2n-i = — i(~ l) n V ; 2(^, z), for broad beam excitation the 
amplitudes ipi and ^2 vary slowly with n, and one can 
thus write ^1,2 ± M) = ipi,2(n,z) ± (dipi^/dn) and 
consider n = £ as a continuous variable rather than as 
an integer index. At the input plane z = 0, the ampli- 
tudes ^i(^,0) and ^2(^,0) are proportional to G(2na) 
and G(2na — a) ~ G(2na), respectively, so that one can 
assume ^i(£,0) ~ ^2(^,0) as an initial condition. Under 
such assumptions, from Eqs.(2) it readily follows that 
the two-component spinor ip(£,z) = (ipi^2) T satisfies 
the one-dimensional Dirac equation 

id z i/j + iaad^ - 5f3^ = 0, (4) 

where 

■=(: ;)."=(* ..)• < 5 » 

Note that a and (3 coincide with the <j x and <j z Pauli 
matrices, respectively. Assuming the normalization con- 
dition j c^d^il 2 + 1^2 1 2 ) = 1, after the formal change 
a — > c, (5 — > mc 2 /h, £ — > x and 2: — > £, Eq.(4) cor- 
responds to the one-dimensional Dirac equation for an 
electron of mass m in absence of external fields [2]. The 
temporal evolution of the spinor wave function ^ for 
the Dirac electron is therefore mapped into the spatial 
evolution of the modal amplitudes ip\ and ^2 in the 
two sublattices A and B. The energy-momentum dis- 
persion relation huj(k) of the Dirac equation (4), ob- 
tained by making the Ansatz ip ~ exp(z/c£ — iuoz) in 
Eq.(4), is composed by the two branches uo(k) = =be(fc), 
corresponding to positive and negative energy states of 
the relativistic free electron, where e(k) = V^ 2 + o~ 2 k 2 . 
Such branches reproduces the two minibands of the bi- 
nary array [see Eq.(3)] near the boundary of the Bril- 
louin zone [see Fig. 1(b)], where k = 2aq — tt. ZB of 
the Dirac electron corresponds to a rapid oscillation 
of the average position (£)(z) = f d£ £(\^i\ 2 + IV^I 2 ) 



around the classical trajectory. The usual method of un- 
derstanding ZB in the framework of the Dirac equa- 
tion is to derive equations of motion for the Heisen- 
berg operators, and show that they oscillate in time 
[1,2]. We instead work directly in the Schrodinger pic- 
ture and consider wave packet evolution in momentum 
space, namely we set ^1,2 (£, z) = J dk$i^(k, z) exp(z/c£) 
and calculate the spectra ^1, 2(^,2) by solving Eq.(4) 
in momentum space. One then obtains ^l^ik^z) = 
G(fc)[cos(ez) =F i(±ak + 5) sin(ez)/e], where G{k) = 
(l/27r) j d^G(2a^) exp(— ik£) is the Fourier spectrum of 
the exciting broad input beam. The average position (£) 
is then calculated as (£)(z) = 2iri j dk[^l{k)dk^i{k) + 
, 0|(/c)9/ c / 02(^)], which yields after some algebra 

<£>(*) = <0(0)+4tt(7 3 2; J ' dk(k/e) 2 \G(k)\ 2 + 

+ 27ra5 2 J dk(l/e 3 )sm(2ez)\G(k)\ 2 (6) 

The last oscillatory term in Eq.(6) is ZB, superimposed 
to the straight trajectory defined by the first two terms 
on the right hand side of Eq.(6). For G(k) spectrally 
narrow at around k = 0, the frequency of ZB is equal 
to 2e(k = 0) = 28; spectral broadening of G(k) is re- 
sponsible for damping of ZB. It should be noted that 
(£), calculated as the average position for the Dirac 
equation (4), basically reproduces the evolution of the 
beam center of mass (n) in real space, defined as (n) = 
En n \ c n\ 2 )/ J2 n \ c n\ 2 • In fact, it is straightforward to 
show that 

(n) = 2(0 + i - 4ir5a J dk(k / e 2 )\G(k)\ 2 sin 2 (e^). (7) 

The last term in Eq.(7) is usually negligible for a spec- 
trum G(k) narrow at around k = (see the examples to 
be discussed below), so that one has (n) ~ 2(£) + 1/2. 
As an example, Fig. 2 shows typical evolutions of field 
intensities |c n (z)| 2 [Figs. 2 (a) and (c)], and correspond- 
ing behavior of beam center of mass (n)(z) [Figs. 2(b) 
and (d), solid curves], for a = 1 and for two values of 
detuning S. In both cases, the array has been excited 
by a broad Gaussian beam G{x) = exp[— (x/wo) 2 ] with 
wo /a = 12, and the results have been obtained by numer- 
ical integration of coupled-mode equations (2) with ini- 
tial conditions c n (0) = i n exp[— (na/wo) 2 ). A clear trem- 
bling motion of the beam, corresponding to ZB, is ob- 
served, with an oscillation frequency (amplitude) which 
increases (decreases) as 5 increases, according to Eq.(6) 
[compare Fig. 2(b) and Fig. 2(d)]. A similar trembling mo- 
tion would be observed for a binary lattice with asym- 
metric intesrite coupling [20]. Damping of ZB oscilla- 
tions, which arises from the broadening of the spectrum 
G(fc), is also visible in Figs. 2(b) and (d); note also that 
the averaged beam path has a nonvanishing drift veloc- 
ity which arises from the second term on the right hand 
side of Eq.(6). 

We finally checked the correctness of the analysis, based 
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Fig. 2. (Color online) (a) Evolution of |c n (z)| 2 in a bi- 
nary array excited by a broad Gaussian beam for a = 1 
and S = 0.6, and (b) corresponding behavior of beam 
trajectory (n)(z) (solid curve) as obtained by numeri- 
cal analysis of Eqs.(2). The dotted curve in (b), almost 
overlapped with the solid curve, reproduces the behavior 
of 2(£)(z) + 1/2, where (£)(z) is calculated according to 
Eq.(5). (c) and (d): same as (a) and (b), but for S = 1. 
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Fig. 3. (Color online) (a) Beam propagation (snapshot 
of \E(x, z)\ 2 ) in a binary array, as obtained by numerical 
analysis of Eq.(l), for broad Gaussian beam excitation 
tilted at the Bragg angle; parameter values are given in 
the text, (b) Corresponding behavior of (x)(z). 

on the tight-binding model (2), and the feasibility of 
an experimental observation of ZB by numerical sim- 
ulations of the wave equation (1) for parameter values 
that typically apply to binary arrays realized in fused 
silica by femtosecond laser writing [17]. A typical nu- 
merical result is shown in Fig. 3 for parameter values 
A = 633 nm, n s = 1.42, a = 10 //m, w\ = W2 = 3.5 /xm, 
Ani = 0.003, An2 = 0.00297 and for an input Gaussian 
beam E(x, 0) of spot size wo = 80 /im, tilted at the Bragg 
angle c± 0.64°. The trembling motion of the beam as it 
propagates along the 10-cm-long array is clearly visible, 
and should be easily observed by microscope fluorescence 
imaging [17]. The beam trajectory in Fig. 3(b) has been 
computed as (x)(z) = j dx x\E\ 2 / j dx\E\ 2 , where the 
integration interval is limited by the size of the numeri- 
cal domain. Note that the first oscillations of (x) internal 
to the shaded area of Fig. 3(b) do not correspond to ZB, 



rather they arise because of an initial beam break up 
and appearance of higher-order beams [as indicated by 
the arrows in Fig. 3 (a)] belonging to higher-order bands 
of the array. Such higher-order beams, however, refract 
at large angles and, after few centimeter propagation, 
they are no more overlapped with the main beam un- 
dergoing ZB. 

In conclusion, a photonic analogue of the trembling mo- 
tion of Dirac electrons has been proposed for spatial 
beam propagation in binary waveguide arrays. As com- 
pared to previous proposals [11, 13], the easy experi- 
mental visualization of beam dynamics in waveguide ar- 
rays [17] should greatly facilitate the way toward the first 
observation of ZB. 
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